Electronic Excitations from a Perturbative ~L~DA-\-GdW Approach 
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We discuss an efficient approach to excited electronic states within ab-initio many-body perturba- 
tion theory (MBPT). Quasiparticle corrections to density-functional theory result from the difference 
between metallic and non-metallic dielectric screening. They are evaluated as a small perturbation 
to the DFT-LDA band structure, rather than fully calculating the self energy and evaluating its 
difference from the exchange-correlation potential. The dielectric screening is desribed by a model, 
which applies to bulk crystals, as well as, to systems of reduced dimension, like molecules, surfaces, 
interfaces, and more. The approach also describes electron-hole interaction. The resulting elec- 
tronic and optical spectra are slightly less accurate but much faster to calculate than a full MBPT 
calculation. We discuss results for bulk silicon and argon, for the Si(lll)-(2x 1) surface, the SiHU 
molecule, an argon-aluminum interface, and liquid argon. 



I. INTRODUCTION 

Many-body perturbation theory (MBPT) has be- 
come the state-of-the-art for excited states in electronic- 
structure theoryji^ Starting from a density-functional 
theory (DFT) calculation, the GW method^ and its 
combination with the Bethe-Salpeter equation^ (BSE) 
allow to investigate the spectra of electrons, holes, 
and correlated electron- hole pairs. The great success 
of MBPT is based on the systematic incorporation of 
Coulomb interaction and polarization effects on all length 
scales, which is not considered in most other electronic- 
structure approaches. The significant computational cost 
of MBPT, however, still constitutes a major obstacle for 
the widespread use of the method. This holds in particu- 
lar for larger-scale systems, like defects, hybrid systems, 
adsorbates, nanostructures, and others. In this paper we 
propose a dramatic reduction of the computational re- 
quirements of MBPT. As a result, the excellent precision 
of standard GW and GW+BSE calculations is slightly 
reduced, but instead the treatment of much larger sys- 
tems becomes possible, thus allowing the investigation of 
spectroscopic features that might be inaccessible other- 
wise. 

As key ingredient we exploit the observation that for 
many systems MBPT, when carried out by (wrongly) 
assuming metallic dielectric screening, approximately re- 
produces the band structure of the underlying DFT cal- 
culation (when employing the local-density approxima- 
tion, LDA). This had already been observed by Wang 
and Pickett^ as well as, by Gygi, Baldereschi, and 
Fiorentinisii^ and was subsequently exploited for model 
QP calculations for various materials/'- As illustration, 
Fig. [1] shows quasiparticle (QP) corrections for silicon 
(Si) and solid argon (Ar). The open circles (o) result from 
a conventional GW calculation (with standard RPA di- 
electric screening, "GVF/RPA"), yielding the well-known 
opening of the band gap (by 0.7 eV for Si and 6.1 eV for 
Ar). The squares (□), on the other hand, come from GW 
calculations which employ metallic screening; these QP 
shifts are close to zero (at least for states near the Fermi 
level). The "metallic" dielectric screening is simulated 
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FIG. 1: (Color online) Quasiparticle corrections of (a) bulk Si 
and (b) bulk Ar. The open circles (o) denote standard GW 
data within RPA screening. The squares (□) result from a 
standard GW calculations, but based on a metallic dielectric 
model function (see text). The filled circles (•) result from 
the present LDA+GdW approach, employing the same basis 
and band-summation details as in the standard GW data. 
The asterisks (*, "fast") result from LDA+GdW with 9 plane 
waves (15 plane waves for Ar) and 8 bands, only (see text). 



by a dielectric model function^— Here we use a model 
based on that of Bechstedt, Enderlein and Wischnewski, 9 
slightly modified (see Sec. Ill B|) for broader applicability. 
Such models are controlled by a few parameters, most 
importantly by the macroscopic dielectric constant, e^. 
Setting eoo=oo turns the screening into that of a metal. 

If the GW method with metallic screening, W meta h 
reproduces the DFT-LDA band structure, one can arrive 
at the true QP band structure by adding a self energy 



AE := iGi(W-W v 
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to the DFT-LDA Hamiltonianir— . This procedure, which 
we label "LDA+GdW" throughout the paper (see be- 
low for details), yields the filled circles (•) in Fig. [T] 
These data agree fairly well with the GW /RP A data, 
at least in the important region near the fundamental 
gap; as a rule of thumb, we find that the gaps from 
AE := iG\(W — W metal) agree with experiment to about 
10 %, which largely corrects the DFT-LDA band-gap 
error of 30-50 %. Most importantly, long-range polar- 
ization effects are included, which is completely missing 
from the short-sighted DFT-LDA. This systematic im- 
provement due to AE is much more important than plain 
agreement of band-structure data with experiment. 

While Eq. (P) describes single-particle states of elec- 
trons and holes, a straight-forward extension to cou- 
pled electron-hole pairs and their optical response is 
easily possible by calculating the electron-hole interac- 
tion (from W) and solving the Bethe-Salpeter equation, 
BSE. 1 ^ Since the self-energy considered in this work is 
in principle still the one resulting from the GW approxi- 
mation (GWA), the corresponding approximations to the 
electron-hole interaction are meaningful in the context of 
LDA+GW, as well. 

As an important consequence of the above findings, 
the calculation of GW-like band structures via AS al- 
lows for a tremendous gain in numerical efficiency in four 
respects: (i) the use of model dielectric functions, (ii) 
small basis-set requirement, (m) small band-summation 
requirement, and (iv) weak influence of dynamical effects. 
The underlying reason for all four issues is that the cur- 
rent approach calculates QP corrections to DFT-LDA as 
a perturbation. Most other GW implementations simply 
replace V xc by J] GW , both of which are in the order of 
magnitude of -10 eV or more. In order to get their dif- 
ference to within 0.1 eV, the GW calculation must be 
carried out with numerical precision better than 1 %. 
Our present approach, on the other hand, starts direcly 
at the QP correction (i.e., AE), which is much smaller 
in magnitude (^1 eV) and much more robust. Here it is 
fully sufficient to evaluate all quantities to within 10 %, 
only, to achieve the same final numerical accuracy in the 
band structure. 

The approach to be proposed in this paper is similar to 
the method by Gygi, Baldereschi and Fiorentini^^ who 
employed their pcrturbative GW method for the com- 
prehensive analysis of bulk crystals. As a key difference 
to their approach, here we employ a different, more gen- 
eral model dielectric function which is flexible enough 
to also describe composite systems containing metals, 
non-metals, molecules, surfaces, interfaces, and more. 
As illustration, we discuss in this paper bulk materials, 
a semiconductor surface, a molecule, a metal-insulator 
junction, and a disordered insulator. In all cases our 
approach yields spectroscopic data of high quality (al- 
though slightly less accurate than the corresponding full 
GW or GW /BSE calculation), demonstrating an appeal- 
ing combination of predictive power, broad applicability, 
and numerical eficiency. 



The paper is organized as follows. In Sec. |TT] we dis- 
cuss the computational approach and the dielectric model 
function employed in this work. In Sec. IIIII characteris- 
tic results for bulk silicon and bulk argon are discussed. 
Sec. IIVI presents results for more complicated systems, 
indicating the potential of the method. The paper is 
concluded by a summary in Sec. |V] 



II. THEORETICAL APPROACH 

In this section we discuss the computational method 
used in this work, including its practical realization and 
underlying physical principles. 

A. Perturbative quasiparticle corrections 

Ab-initio quasiparticle (QP) band structures result 
from the electron self-energy operator E(£'). The 
state-of-the-art approach to E is given by Hedin's GW 
approximation, 3 which is usually evaluated and employed 
on top of an underlying density-functional theory (DFT) 
calculation. The typical procedure employs DFT data to 
generate the single-particle Green function G\ and the 
screened interaction W (usually within the random-phase 
approximation). Thereafter, the resulting self-energy op- 
erator E = iG\W replaces the DFT exchange-correlation 
potential, V xc , arriving at a QP Hamiltonian of 

HQP ~H DFT + iG 1 W-V xc (2) 

This procedure is commonly labelled "many-body per- 
turbation theory" (MBPT). However, this does not mean 
that Eq. (J5|) would be evaluated truly perturbatively in 
the sense that the smallness of the difference {iG\W — 
V xc ) would be exploited. Instead, both terms, iG\W 
and V xc , are evaluated separately, taking their difference 
afterwards. The QP corrections are thus obtained as 
(small) differences between two (rather large) quantities, 
both of which have to be evaluated independently and 
with high precision. Simply speaking, in order to get 
their difference (often ~1 eV) to within 0.1 eV (which 
is the accuracy expected from MBPT), both iG\W and 
V xc (being of the order of ~ —10 eV or more) need to be 
evaluated with a precision of 1%. The underlying reason 
for this problem is the quite different conceptual origin 
of the two terms, which makes it difficult to formulate 
their difference in analytic terms 

Fortunately, there does exist some pragmatic link be- 
tween E and V xc : The self energy of the homogeneous 
electron gas (as a function of the energy of a given state) 
is nearly constant (see Ref. Q) and thus nearly coincides 
with V xc . [Note that this might not be truly fulfilled by 
approximations to E, like the GWA, which might suffer 
from offsets.] This behavior is reflected by the observa- 
tion that in bulk metals, QP corrections (from GWA) to 
DFT-LDA band structures are very smalL 15 ' 16 In other 
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words, DFT (at least within the local-density approxima- 
tion, LDA) does contain correct spectral properties of the 
quasiparticles, at least for homogeneous systems. Within 
LDA, however, these spectral properties are by construc- 
tion still those of the metallic system (jellium) from which 
the LDA exchange-correlation data originate, and this 
metallic behavior (in particular, metallic screening) is a 
built-in property of the exchange-correlation potential, 
even when applied to non-metallic systems. A generaliza- 
tion of this statement would imply that V xc s» iGiW meta i 
(provided that iG\W is a good approximation to S) 
with the appropriate W me t a i (i-e. metallic screening). In 
fact, GW studies employing metallic screening (including 
ours, see Fig. [1]) confirm that V^ c |nk) sa iG iW me tai\nk) 
for most electronic states |nk). 

Based on the working hypothesis that for non- 
homogeneous, non-metallic systems the largest differ- 
ence to metallic behavior is the difference in screening, 
and employing V xc f=s iG\W me tai, one arrives at the QP 
Hamiltonian 



jjQP „ fjDFT-LDA + _ 



(3) 



in which AE = iG\{W — W meta i) acts as a self energy, 
yielding QP corrections (cf. Eq. (JTJ). The most im- 
portant change to Eq. @ is the fact that Eq. no 
longer evaluates the difference between the self energies 
(given by iGiW and V xc ), but the difference in screen- 
ing: (W — Wmetai)- This difference is much simpler and 
faster to treat than the difference between self energies. 
Note that on the other hand, the final accuracy of the QP 
band structure might be less than 0.1 eV because the en- 
tire approach, although being much more efficient from 
a numerical point of view, is based on the assumption 
that V xc ~ iGiWmetah meaning a further approximation 
in addition to the GW approximation. The assumption 
that V xc s» iGiWmetai should be checked carefully for 
each system class. 

As discussed below, the use of Eq. ^ allows for sev- 
eral numerical simplifications (see Sec. Ill D|) . leading to 
a higher efficiency than conventional GW calculations, 
allowing to tackle more complex systems. One of the 
most important facilitations is the use of model dielec- 
tric functions (see next section) instead of calculating the 
screening within the random-phase approximation. 



B. Model dielectric function 

The calculation of the dielectric function within the 
random-phase approximation (RPA), which is the com- 
mon procedure within MBPT, is one of the bottlenecks 
of the method. A simplified evaluation of the dielectric 
function is an important contribution to improving the 
efficiency of MBPT (even without the considerations of 
the previous section). For this reason model dielectric 
functions are sometimes employed to avoid the RPA»i£~— 
In the present context, the use of models is also manda- 
tory for another reason: The key ingredient of the present 



theory is the difference between the correct screening of 
the (non-metallic) system and its (hypothetical) metal- 
lic counterpart. This is only useful and well-defined if 
both types of screening result from the same approach, 
allowing to tune the screening from "correct" to "metal- 
lic" in a seamless manner. It is, however, unclear how 
the RPA could be used to simulate metallic behavior of a 
non-metallic system. An appropriate model is therefore 
a necessity of the current approach. 

Examples are the models proposed by Bechstedt, En- 
derlein, Wischnewski, and Falter and by Levine, Hybert- 
sen and Louie^r— We have tested these models in the 
present context and find that they yield essentially the 
same results as the ones to be discussed below. In their 
original form, however, these models have one signifi- 
cant disadvantage which may hinder their application to 
more complex systems: they were formulated for systems 
that are characterized by one common dielectric constant 
without spatial variation. This makes it difficult to em- 
ploy them for systems in which the screening shows spa- 
tial variation, like interfaces, molecules, etc. 

Instead we propose a model that is based on a combina- 
tion of localized and delocalized quantities. The system 
may consist of N atoms (at positions Tj) in a (periodi- 
cally repeated) unit cell or supercell of volume V, with 
reciprocal lattice vectors G. To each atom we attribute 
a static charge-density response x ( see below) and an 
effective volume Vj . The dielectric function of the whole 
system is then obtained as 



|q + G||q- 



G'| f-f V 

3 = 1 



The volume attributed to each atom controls the weight 
which the atom contributes to the response. The trans- 
formation from the charge-density response to the di- 
electric function further involves a convolution with the 
Coulomb interaction, i.e. the multiplication by l/(|q + 
G||q + G'|) in Eq. Note that we work with a sym- 
metrised dielectric function ! 23 ' 24 

It was suggested by Bechstedt et al. to describe the 
charge-density response of a (homogeneous) system with 
dielectric constant by^ 



f 



Q 2 



1 qfr F {p) ujp(p) 



(5) 



where the Thomas-Fermi wave number qtf and plasma 
frequency uip depend on the average electron density p. 
Eq. ((SJ) is related to the Lindhard dielectric function. In 
combination with Eq. J1J (for homogeneous systems, dis- 
regarding the summation over atoms), x = Q 2 '.f would de- 
scribe the dielectric function. In particular, for Q— >0 one 
would correctly obtain e(Q)— >e oa (if e QO <oo, i.e. for non- 
metals) or e(Q) — > 1 + Qtf/Q 2 e oo=oo, i.e. for metals), 
respectively. The generalization to non-homogeneous 
systems is less clear. While the large- Q behavior, which 
reflects the short-range reaction of electronic charge to 
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external fields on the sub-atomic length scale, appears 
realistic for non-homogeneous systems as well, the real- 
istic incorporation of atomic-length-scale charge-density 
variation and of local fields (i.e. off-diagonal matrix ele- 
ments of Xgg'W) ^ s l ess c l ear - 

Here we propose to model the charge-density response 
attributed to each atom by: 

atom with metallic response: 
xg G <(q) = V //(|q+G|;^,oo)/(|q + G'|;p J -,oo)- 

•|q + G||q + G'| • e -7 3 (G'-G) 2 el (G'-G)r J (g) 
atom with nonmetallic response: 
xg G '(q) = ^/(Iq+Ghp^e.O/dq + G'l;^^,). 

•(q + G)(q + G') • e -7 3 -(G'-G) 2 ei (G'-G)T, (?) 

In both cases, the factor \//(|q + G|)/(|q+ G'|) is a 
reasonable average of the large-Q behavior in directions 
(q + G) and (q + G'). The phase factor for each atom 
results from the position of the atom within the unit cell 
or supercell. The factor exp[— jj(G' — G) 2 ] describes the 
spatial extent of the charge density of atom j. With- 
out this factor (or with jj— >0), the model describes a 
sharp point-charge-density response at position Tj . With 
7j— >oo all local fields would be switched off, turning 
the model into that of a homogeneous system again. 
For a non-zero, finite value of 7^, the factor exp[z(G' — 
G)Tj] exp[— 7j ;(G' — G) 2 ] is the Fourier transform of a 
Gaussian-shaped charge density ~ exp[— (r — Tj) 2 / (47,)] 
centered at Tj. In short, this means that the charge- 
density response is neither perfectly local (i.e., exactly 
at Tj) nor completely delocalized (except for a truly ho- 
mogeneously system, to be characterized by jj — > 00). 
Instead, the charge-density response of an atom origi- 
nates from its charge density (or at least from that of 
the polarizable electronic states), and its spatial form is 
included in the model. Correspondingly, approxi- 
mates the radius of the atom. It should be noted that the 
term exp[— 7j(G' — G) 2 ] • exp[z(G' — G)tj] corresponds 
to the factor p(G — G')/p(0) (i.e., Fourier transform of 
the charge density) in the model by Bechstedt et al£ 
In our model the charge density of the entire system is 
approximated by a composition of atomic contributions 
with simplified shape. 

A particular role is played by the factors |q+G||q+G'| 
(for metallic response) and (q + G) • (q + G') (for non- 
metallic response) in Eqs. (JSJ) and ([7]). These factors 
reflect the qualitatively different origin of the response 
of metallic and non-metallic systems. For metals, long- 
range charge fluctuations and displacements are possi- 
ble, resulting from intraband transitions near the Fermi 
level. Such displacements lead to charge accumulation 
at some atoms and charge depletion at others. Here 
our model assumes that such charge accumulation or de- 
pletion would basically show the same spatial structure 
as the original charge density of the atom (modeled by 
exp[-(r - -r J ) 2 /(4 7i )]), i.e. 5 Pj (r) ~ Pj (r). 



The charge-density response of a non-metal, on the 
other hand, is of completely different origin. Here the 
response to an external field is mainly given by a short- 
range displacement of charge density from one side of 
the atom to the other, i.e. by a polarization of the 
atom. In many cases, this polarizability is dominated 
by transitions from s orbitals to p orbitals or vice versa. 
The spatial structure of such si-tp polarizability is given 
by a factor (r— Tj) ■ (r'—Tj), leading to a factor of 
(q+G)-(q+G') in reciprocal space. Again, the additional 
factor exp[— (r— Tj) 2 /(Ajj)] (or its reciprocal-space coun- 
terpart) reflects the fact that the response comes from the 
whole atom (including some spatial extent) rather than 
from a single point. 

The model can also be generalized to the case of 
anisotropic response, e.g. if an atom is embedded in 
a non-isotropic chemical environment, like in molecules, 
at surfaces, in sp 2 -bonded carbon, in atomic monolayers 
on a substrate, or similar. The same holds for mate- 
rials with an anisotropic dielectric-constant tensor. In 
both cases, the response of each atom should be modeled 
with a direction-dependent dielectric-constant parameter 
tj (q) . Since such a situation can be expressed in terms of 
the three principal axes n^-* and corresponding principal 
values e^ h > of a 3x3 tensor (see next section for details), 
a straight-forward generalization of Eq. (J7J is possible: 

atom with nonmetallic response: 

xg G '(q) = 

^^//(Iq+Gh^-.ef )/(|q + G'|;p J , £ f ) )• 

k=l 

•[(q + G)nf ] • [nf (q + G')] . e -^(G'-G)y ( G'-G^ 

One especially useful feature of the model proposed 
in this work is the possibility to combine metallic and 
non-metallic response in one system. This is particularly 
relevant for adsorbates on metallic substrates, for metal- 
insulator interfaces etc. Here our model simply allows 
to attribute metallic parameters (i.e., Eq. ©) to some 
atoms and non-metallic parameters (i.e., Eqs. ([TJ or ©) 
to others. For the construction of W me t a i, finally, we 
simply take metallic response of all atoms (i.e., Eq. ©). 
Apparently, for metals (or metallic regions) the dielectric 
function is the same in both cases. 



C. Determination of the model parameters 

The determination of the parameters is a particular 
task. Fortunately, the final use of the model for differ- 
ences between metallic and non-metallic screening makes 
the entire approach insensitive to the actual choice of the 
parameters Vj, pj, and jj. Within this work, we simply 
attribute a realistic volume to each atom (for silicon, e.g., 
we choose Vj — 20 A 3 , which is the volume per atom of 
bulk Si), as well as a realistic valence electron number 
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(for silicon, apparently Nj—4). The average electron den- 
sity pj = Nj/Vj defines the Thomas- Fermi wave number 
Qtfj and plasma frequency wpj- for this atom, to be used 
in f(Q). The parameter 7^ is obtained from least-square 
fitting of the atomic charge density by a Gaussian func- 
tion. These parameters are used for both the metallic 
and the non- metallic response of atom j. 

For non-metals one needs the dielectric-constant pa- 
rameter ej (or the principal axes and values of the corre- 
sponding tensor for anisotropic situations). Such values 
can either be taken from experiment, or they are calcu- 
lated for the particular system. One possibility is given 
by the evaluation of the small-q limit of £G=o,G'=o(q) 
from the electrical-dipole operator applied to the inter- 
band transitions of the system, leading to a 3x3 tensor 
of 



Energy cutoff [Ryd] 
2 



occ empty 
v c k 



'vk\p a \ck){vk\p b \ck) 



(£ C k 



(9) 



(a,b = x,y,z) from which the principal axes and values 
can be evaluated. Note that Eq. © does not contain 
local-field effects. However, since the resulting ej enter 
our model before the inversion of £g,g'(i) (which then 
leads to the local- field effects) , the employment of local- 
held-free parameters is not a problem, but rather a re- 
quirement of the model. 

In many systems the responses of the various atoms 
will differ from each other, leading to the question of 
distributing the results of Eq. (|9]) over the individual 
atoms. Here we propose to employ an atom-centered 
local-orbital basis for the calculation of the electronic 
states, \rik). Such a basis allows to decompose the dipole 
matrix elements (wk|p a |ck) into individual contributions 
of each atom j, i.e. one can focus on atom j and switch 
off the dipole strength of all other atoms. In this case 
Eq. ([9]) yields an individual result for each atom alone, 
allowing to find out the individual parameters of each 
atom in the system. 



D. Numerical efficiency 

As an important consequence of the above findings, 
the calculation of GW-like band structures via AS al- 
lows for a tremendous gain in numerical efficiency in four 
respects: (i) the use of model dielectric functions, (ii) 
small basis-set requirement, (Hi) small band-summation 
requirement, and (iv) weak influence of dynamical effects. 
As mentioned, the underlying reason for all four issues is 
that the current approach calculates QP corrections to 
DFT-LDA as a perturbation. The four issues of efficiency 
(i)-(iv) deserve detailed discussion. 

(i) The advantage of working with a dielectric model 
function rather than employing the random-phase ap- 
proximation has already been pointed out in the last sec- 
tion. In particular, the perturbative idea of Eqs. (|T|) and 
requires a model, even beyond the issue of numerical 
efficiency. 
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FIG. 2: (Color online) Gap energies of Si (indirect mini- 
mum gap, direct gap at L, and direct gap at X), calculated 
within a full GW calculation (employing RPA screening) and 
within the present LDA+GriVF scheme, (a) Dependence of 
the gap energies on the plane-wave basis representation of W 
(or (W — Wmetai), respectively), as controlled by the energy 
cutoff (upper axis), (b) Dependence of the gap energies on 
the number of bands considered in the band summation in- 
side the self-energy operator. [ Note that this does not refer 
to the number of bands considered in the calculation of the 
RPA screening.] 



(ii) The basis-set requirements for (W — W meta i) are 
much weaker than for W alone for two reasons: (ii.l) 
Within full GW the bare-exchange contribution requires 
a large basis for convergence. In our present approach, 
on the other hand, the bare-exchange effects are the same 
in GW and GW meta i and thus cancel each other, (ii.2) 
While both W and W meta i are structured in real space, 
their difference is a rather smooth function and converges 
with very few basis functions. Fig. [5] a shows repre- 
sentative gap energies of Si as a function of the plane- 
wave basis size used for W. Full GW requires about 
60 plane waves (~5 Ryd cutoff) for reasonable accuracy 
while the LDA+GdW data are already converged with 
9 plane waves (~ 1.4 Ryd). Similarly, basis-set conver- 
gence for Ar requires cutoff energies of about 15 Ryd for 
GWA, but only about 2 Ryd for LDA+GdW . 

(Hi) The band summation in AS is less demanding 
than in a full GW calculation because the influence of the 
higher conduction bands (via Gi) is weak. This is shown 
in Fig. [2]b. The use of about as many conduction bands 
as valence bands is sufficient for LDA+GdW (at least for 
states near the gap), while conventional GW calculations 
are usually performed with at least about 10 times more 
conduction than valence bands. This behavior again re- 
sults from the smooth spatial structure of (W — W me tai)- 
To summarize statements (ii) and (Hi), Fig. [1] includes 
data from 1.4 Ryd cutoff (2.0 Ryd for Ar) and four con- 
duction bands in G\ as asterisks (*, "fast"). The agree- 
ment with the converged LDA+GdW data shown by (•) 
is sufficient, except for higher-energy states. 

(iv) Within conventional GW calculations, the corre- 
lation part of S (which is subject to dynamical effects) 
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can be as large as 5-10 eV. Our AS, on the other hand, is 
much smaller in magnitude and thus much less sensitive 
to dynamical effects. This allows to treat AS on the level 
of the static COHSEX approximation^ which we em- 
ploy in all LDA+GdW calculations in this paper. Note 
that this does not apply to the GW /RP A and GW /Metal 
reference calculations in this paper, all of which include 
dynamical effects by using a plasmon-pole approxima- 
tion. In those cases, the generalization of the static model 
of Sec. Ill Bl to a dynamic dielectric function is realized 
by evaluating the /-sum rule^ 



III. RESULTS FOR BULK SILICON AND 
ARGON 

The QP corrections to DFT-LDA for bulk Si and 
Ar ar compiled in Fig. [TJ As discussed above, full 
GW calculations using a metallic W yield QP correc- 
tions close to zero, opening the possibility of perturba- 
tive LDA+GdW as proposed in Sec. [TTJ In fact, the 
LDA+GdW data are close to those of a full GW cal- 
culation employing correct, non-metallic screening from 
RPA (open circles). There are, however, some deviations 
(related to the LDA+GdW method as such, and also 
to the dielectric model function). For Si, for example, 
the lowest valence bands observes very small QP cor- 
rections within GW /RPA, but significant negative QP 
corrections within LDA+GdW . Furthermore, the QP 
corrections for the conduction bands appear to be less 
accurately reproduced by LDA+GdW than for the up- 
per valence bands. Additional deviations are observed 
for the "fast" LDA+GdW approach (at minimal basis-set 
and band-summation specification), in particular for the 
higher conduction bands. These details nonwithstand- 
ing, we can conclude that LDA+GdW yields sufficient 
accuracy if one is interested in states near the funda- 
mental gap. Furthermore, systematic deviations between 
LDA+GdW and GW /RPA can be expected to be simi- 
lar in bulk systems and other, more complicated systems 
of the same material (like, e.g., surfaces, nanostructured 
systems, interfaces, etc.). LDA+GdW will allow for sys- 
tematic comparison between the spectral data of such 
systems. 

The quite reliable LDA+GdW band structures can be 
employed to yield reasonable optical spectra, as well. To 
this end we include electron-hole interaction on the level 
of the Bethe-Salpeter equation^ The interaction kernel is 
calculated with the same (non-metallic) dielectric model 
function and same basis as the band structure. One ex- 
ception is the unscreened exchange interaction between 
electrons and holes (originating from the Hartree poten- 
tial) which may require a larger energy cutoff than the 
screened interaction and must be treated separately. Its 
calculation is relatively cheap and does not affect the ef- 
ficiency of our approach. 

Fig. |3] shows the macroscopic imaginary dielectric 
function £2(w) of bulk Si and Ar. The thin dotted curves 
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FIG. 3: (Color online) Optical spectrum (imaginary part of 
macroscopic dielectric function) of bulk Si and of bulk Ar. 
Experimental data are from Ref. [27l - [29l (Si) and from Ref. [3(] 
(Ar). Note that spin-orbit coupling (leading to the measured 
double-peak structure of the Ar exciton) is not included in 
our calculations. 



display the spectrum from the uncorrelated LDA in- 
terband transitions, which is usually qualitatively and 
quantitatively wrong (in particular for insulators, like 
Ar). The dashed lines are reference data from a full 
GVF+BSE/RPA calculation, which can be considered as 
the state-of-the-art approach to e2(w). The solid lines 
display our current results, including the drastic numer- 
ical simplification ("fast") as outlined above. In com- 
parison with the GVF+BSE/RPA data and with exper- 
iment, the LDA+GdW data are very gratifying. They 
correctly yield the two characteristic peaks (at 3-3.5 eV 
and at 4-4.5 eV) of the Si spectrum. For Ar, we obtain 
an exciton peak at 12.2 eV (GIF+BSE/RPA) or 12.1 
eV (LDA+GdW), respectively. In experiment, the spin- 
orbit interaction (neglected in our present work) splits 
the exciton into two peaks at 12.0 and 12.2 eV. Fur- 
thermore, a second excitonic peak is found near 13.5 
eV. The slight deviations between LDA+GdW and the 
GPF+BSE/RPA reference data mostly result from corre- 
sponding deviations in the band structure. Note that for 
Ar, the agreement is even better than can be expected 
from our current approach. Most importantly, the differ- 
ences between LDA+GdW and GW+BSE/KPA arc not 
significantly larger than the deviations from experiment, 
thus advertising LDA+GdW as a useful alternative. 
Compared to the LDA interband spectrum, a tremen- 
dous improvement of explanatory power is achieved. 
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[eV] LDA GW/ GW / LDA+ Exp. 

Metal RPA GdW 

bulk E^ p n 0.49 0.51 1.10 0.95 1.17 a 

D up (J) 0.0 0.0 0.0 0.0 0.0 6 

DdcwnjJ) 0.4 0.4 0.7 0.8 0.7 C 

a Ref. H3 6 Ref. |H c Ref. M 

TABLE I: Characteristic band-structure data for the Si(lll)- 
(2x1) surface, which is dominated by two dangling-bond 
states derived from the Pandey-chain termination. At the 
J point of the surface Brillouin zone the related bands (oc- 
cupied D up and unoccupied Ddown state) are closest to each 
other and define the surface band gap. 34 



[eV] 
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GW 1 
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Exp. 
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[35J 


Elomo 


-13.5 
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-17.8 
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-8.4 


-9.0 


-12.5 


-11.8 


-12.6 


Elumo 


-0.6 


-0.2 


0.4 
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^triplet 






8.0 


7.6 




^singlet 






9.0 


8.3 


8.8 



TABLE II: Spectral data of the S1H4 molecule, which is dom- 
inated by quantum confinement and shows the typical elec- 
tronic excitations of a small molecule. 



IV. RESULTS FOR MORE COMPLEX 
SYSTEMS 

We have tested the LDA+GdW approach for a num- 
ber of inhomogeneous systems, starting from the bulk 
materials (Si and Ar) discussed above. 

A. Si(lll)-(2xl) surface and silane molecule 

Based on the experience with bulk silicon, we investi- 
gate two prototypical systems of silicon in reduced dimen- 
sions, i.e. the SiH.4 molecule and the Pandey-chain ter- 
minated Si(lll)-(2xl) surface. Both systems have been 
intensively studied in theory and experiment (see, e.g., 
Ref. H3 and [3^ and references therein). Here we focus on 
their electronic structure within the present LDA+GdW 
approach. 

In the case of the Si(lll)-(2xl) we focus on the 
band structure of the Pandey-chain derived dangling- 
bond states i 32 ' 33 i 37 ~ — The dangling bonds result from 
the lower coordination (three-fold instead of four-fold) 
of the Pandey-chain atoms, leading to one occupied and 
one empty state within the bulk band gap.— These two 
bands constitute one of the most intensively studied sur- 
face electronic structures. At the J point of the sur- 
face Brillouin zone the two bands are well separated from 
the silicon bulk states and define the surface band gap. 
Within LDA, this gap suffers from the same type of band- 
gap underestimation as all semiconductor systems. Here 
we observe a value of 0.4 eV, much smaller than the ex- 
perimental result of 0.7 eV from a combination of direct 
and inverse photoemission (see Tab. [I])] 32 ' 33 



Within GW /RPA, the surface bands are significantly 
shifted and result in very good agreement with the mea- 
sured data.— ~ 34 ' 40 i 41 It is most gratifying to see that 
this behavior is also given by the present LDA+GdW-^ 
approach, which yields a surface gap energy of 0.8 eV. 
This good agreement also holds for the absolute ener- 
getic position (with respect to the bulk band structure). 
Both for the occupied and for the empty band, the data 
from GW /RPA, LDA+GdW, and experiment all agree 
to within 0.1 eV. 

The screening properties for this calculation have been 
obtained from the approach as outlined in Sec. Ill CI 
yielding individual screening properties for each atom. 
Here we find that the charge-density response of the sur- 
face atoms is slightly larger than that of the bulk-like 
atoms in the center of the slab. The response of the 
bulk-like atoms agrees with that of a true bulk calcula- 
tion to within 10 percent. At the surface, on the other 
hand, the smaller surface band gap, the 7r-conjugated na- 
ture of the Pandey chain, and the vicinity of the vacuum 
lead to an anisotropic response. Perpendicular to the sur- 
face, the response is reduced by about 25 % (leading to a 
dielectric-constant parameter of about e^- =9 instead of 
the bulk value of ej=12). Parallel to the Pandey chain, 

on the other hand, the response is doubled to e^ ll) =24. 

As another, even more extreme case for silicon 
in reduced dimension, we discuss the silane molecule 
(SiH4) j 35 ' 36 i 43 Its electronic structure is dominated by 
quantum confinement, leading to much larger band gaps 
and QP corrections than for extended semiconductors. 
All relevant data are compiled in Tab. |TTJ Compared 
to the LDA data, the occupied states (i.e. the lowest 
occupied molecular orbital, LOMO, and (three-fold de- 
generate) highest occupied molecular orbital, HOMO), 
are lowered in energy by more than 4 eV. Here the cur- 
rent LDA+GdW approach reproduces these QP shifts to 
within about 0.5 eV. The lowest unoccupied molecular 
orbital (LUMO), on the other hand, is shifted to higher 
energies by 1.0 eV (GW/RPA) or 1.2 eV (LDA+GdW), 
respectively. Based on these reliable data for single- 
particle states, LDA+GdW also yields reasonable data 
for charge- neutral electron-hole excitations (see Tab. ILT1) . 
Here we take the lowest-energy singlet and triplet ex- 
citation as representative examples. While GVK+BSE 
within RPA yields data in excellent agreement with 
experiment 3 ^ and with other many-body and quantum- 
chemical methods^ the data from LDA+GdW show 
slightly lower excitation energies. The deviations are in 
the order of 0.5 eV and correspond to the differences in 
the band-structure energy of the HOMO state, for which 
LDA+GdW yields a slightly too high value. Neverthe- 
less, in light of the huge QP corrections and very strong 
electron- hole interaction of about 5 eV in SiH4, we con- 
sider the accuracy of LDA+GdW (i.e. yielding QP shifts 
and electron-hole binding to within 20 %) extremely grat- 
ifying. 

Similar to the case of the Si(lll)-(2x 1) surface, screen- 
ing in SiH4 differs significantly from that of bulk silicon. 
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FIG. 4: (Color online) Local density of states in a mono- 
layer of Ar, combined with a 5-layer aluminium(OOl) slab in 
a heterostructure. 

The much larger gap reduces the charge-density response 
strongly. Our approach of locally evaluating the density- 
response contribution of each atom yields an isotropic 
response of the silicon atom to be described by €^=3.75 
(and similar results for the H atoms), i.e. weaker than 
bulk Si by a factor of 4. Such strong reduction for chem- 
ically passivated silicon in confined geometries was al- 
ready found earlier 4^ 

We close this section by mentioning that for both sys- 
tems, Si(lll)-(2x 1) and S1H4, the underlying reason for 
the success of the LDA+GdW approach is again given by 
the reproduction of the DFT-LDA band-structure data 
when metallic screening is employed in a full GW calcula- 
tion. The corresponding data are included in Tab. [I] and 
Tab. [TTJ In particular for the Si(lll)-(2x 1) surface, this 
mandatory condition for the applicability of LDA+GdW 
is nearly exactly fulfilled. For the SiH4 molecule some 
difference of the order of 0.5 eV are found. Considering 
the massive deviation of this system from a homogeneous 
metal, even this agreement to within 0.5 eV is an amazing 
result. 



B. Argon systems 

Spatially varying dielectric response is also present in 
metal-insulator heterostructures. As an example Fig. 
U shows the single-particle spectrum of a periodic het- 
erostructure composed of five atomic layers (10 A) of 
aluminium and one atomic layer (3 A) of argon, stacked 
along the Al(001) direction. This system combines metal- 
lic screening in Al with insulating behavior in Ar, which 
has significant consequences on the QP energetics 
Here we focus on the local density of states (LDOS) in 
the Ar monolayer. The LDOS between -12 eV and -6 eV 
results from the upper valence states of Ar (3p) , while the 
LDOS above 2 eV comes from the Ar conduction bands, 
with increasing admixture of Al states at higher energy. 
The LDOS inside the Ar band gap (-6 eV to +2 eV) re- 
sults from spill-out of Al states into the Ar layer. The 
most interesting feature is the rather small QP correc- 
tion of the argon states, which (in GW /RPA) amounts 
to -1.7 eV (+0.2 eV) for the upper valence (lower con- 
duction) states, yielding a total correction of 1.9 eV for 
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FIG. 5: (Color online) Optical spectrum of the excitons in 
non-crystalline argon (from molecular dynamics of 64 atoms 
at 300 Kelvin, at the solid-argon density — see text), result- 
ing from the current LDA+GdW approach. The dashed line 
indicates the position of the exciton in the periodic crystal. 
The inset shows the Ar-Ar pair-correlation function of the 
MD simulation. 



the fundamental gap of Ar. In bulk Ar, on the other 
hand, the gap-edge states observe QP shifts of -4.1 eV 
and +2.0 eV, yielding a gap correction of 6.1 eV (cf. Fig. 
[TJ . The presence of metallic screening in the immediate 
neighborhood significantly weakens the QP shifts due to 
image-state effects^ - — both for holes and for electrons 
(by about 2 eV each). It is most gratifying to see that in 
our LDA+GdW approach (again with a plane-wave cut- 
off of 2 Ryd), the spatial set-up of the dielectric model 
function (cf. Eq. (U)) reproduces these effects. Here 
LDA+GdW yields QP shifts of -1.8 eV and +0.2 eV for 
the band-edge states and a gap correction of +2.0 eV, 
compared to the Ar bulk data (see Fig. [1} of -3.5 eV, 
+2.8 eV, and +6.3 eV. We conclude that the LDA+GdW 
approach is a suitable method for addressing electronic 
properties of metal-nonmetal junctions. 

As a last example for the potential of our method, 
Fig. [5] shows the exciton spectrum of non-crystalline ar- 
gon. At zero temperature argon forms a periodic face- 
centered cubic (fee) lattice, which can easily be treated 
by MBPT (see Sec. [II]), leading to the results as discussed 
in Sec. IIII1 For this periodic solid the exciton yields a 
sharp line (except for dynamical broadening effects from 
self trapping or similar, that are completely neglected 
here). This changes in the case of non-periodic argon, 
like in its liquid or amorphous state. Such systems may 
be described by sufficiently large supercells. At present 
we investigate the spectra resulting from a 64-atom cell 
(consisting of 4x4x4 fee unit cells) and exploit its spec- 
tral features from the T point of the supercell, only. For 
the periodic fee crystal this yields an exciton at 12.01 eV 
excitation energy (slightly lower than the value reported 
in Sec. IIII1 which was obtained from the standard fee 
unit cell containing one atom, and 500 k-points in the 
BSE). Within this configuration (which is computation- 
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ally much more demanding than a simple one-atom-fcc 
calculation and is extremely demanding for the standard 
GW +BSE Hamiltonian) the spectrum of liquids or amor- 
phous systems can be evaluated. At present we simply 
consider argon at its solid-state density (for compari- 
son sake), but heated to 300 Kelvin (although this is an 
unrealistic high temperature for argon at this density). 
We perform a constant-temperature molecular-dynamics 
simulation (using a simple Lennard- Jones interatomic po- 
tential), leading to the Ar-Ar pair-correlation function 
shown in the inset of Fig. [5] Such a simulation is cer- 
tainly not fully realistic in terms of describing liquid or 
amorphous systems; nonetheless it yields structural el- 
ements that may very well be present in liquids. The 
pair-correlation function clearly exhibits structures be- 
yond harmonic vibrations (like, e.g. the vanishing of the 
second-nearest-neighbor peak at 5.3 A), thus prohibiting 
a perturbative electron-phonon interaction treatment in 
the evaluation of the spectrum. Instead, our LDA+GdW 
approach (averaged over 20 snapshots of the MD run) 
easily allows to evaluate the spectrum. The data shown 
in Fig. [5] clearly demonstrate three important features: 
(i) the exciton line is significantly broadened, (ii) the 
broadening is asymmetric, leading to substantial non- 
zero amplitude well above the exciton energy, and (iii) the 
maximum of the peak is at lower energy than in the peri- 
odic system. The third feature is related to the fact that 



in the pair-correlation function, the first maximum also 
occurs at smaller distance (3.4 A) than the fee nearest- 
neighbor-distance (3.7 A), which is a consequence of the 
anharmonicity of the Ar-Ar interatomic potential. 



V. SUMMARY 

In summary, we have discussed an extremely efficient 
modification of standard many-body perturbation the- 
ory (GW method plus Bethe-Salpeter equation). Based 
on the observation that metallic screening in the GW 
method approximately reproduces the DFT-LDA band 
structure (which should be checked for each material), 
quasiparticle (QP) corrections to DFT-LDA are obtained 
in a truly perturbative approach at minimal cost, pro- 
vided that the dielectric screening can be described by 
an appropriate model. The resulting QP band structures 
and optical spectra (including electron-hole interaction) 
are slightly less accurate than those from conventional 
GW+BSE, but they include all Coulomb-interaction ef- 
fects (like screening, electron-hole binding etc.) in a phys- 
ically correct way, allowing to systematically investigate 
excitations beyond DFT and beyond the computational 
limits of conventional MBPT. 
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